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We present a multiple time step algorithm for hybrid path integral Monte Carlo simulations 
involving rigid linear rotors. We show how to calculate the quantum torques needed in the simulation 
from the rotational density matrix, for which we develop an approximate expression suitable in the 
case of heteronuclear molecules. 

We use this method to study the effect of rotational quantization on the quantum sieving prop- 
erties of carbon nanotubes, with particular emphasis to the para-T^I para-Yiz selectivity at 20 K. 
We show how to treat classically only some of the degrees of freedom of the hydrogen molecule and 
we find that in the limit of zero pressure the quantized nature of the rotational degrees of freedom 
greatly influence the selectivity, especially in the case of the (3,6) nanotube, which is the narrowest 
tube that we have studied. 

We also use path integral Monte Carlo simulations to calculate adsorption isotherms of different 
hydrogen isotopes in the interior of carbon nanotubes and the corresponding selectivity at finite 
pressures. It is found that the selectivity increases with respect to the zero pressure value and tends 
to a constant value at saturation. We use a simplified effective harmonic oscillator model to discuss 
the origin of this behavior. 
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I. INTRODUCTION 

Carbon nanotubes have received a lot of attention in the past years in many fields of science and engineering, 
ranging from metallic properties, to adsorption of light gases. 

The availability of these potentially narrow (subnanometer) channels has raised the interest about the possibility of 
using them as isotope sieves. The presence of quantum effects on adsorption has been known for quite a long time 1 ^ 
and Krylov et al£ were the first to show, in the framework of a simple rigid wall model, the possibility of separating 
helium and hydrogen isotopes by using the different zero-point energy of the two confined species. 

The idea was developed by Johnson and collaborator s 4 ! 5 ' 6 , with particular emphasis about the calculation of the 
selectivity of T2/H2 mixtures in the interior channels of carbon nanotubes. These authors used an interaction potential 
model which neglects the molecular structure of hydrogen and developed a method - known as the simple theory - 
to calculate the selectivity in the limit of zero pressure (i.e. neglecting hydrogen-hydrogen interactions) by using the 
single particle energy levels in the confined systemi 4 ^ 5 - 

They conclude that the (3,6) nanotube would be able to show a selectivity of the order of 10 5 at zero pressure and 
20 K and calculated the expected selectivities for a wide set of nanotubes, thus displaying the dependence of this 
quantity on the radius of the tubes. They later extended their calculations, using path integral Monte Carlo methods, 
to finite pressures £ 

Hathorn et alZ- were the first to address the effect of the rotational degrees of freedom of hydrogen on the selectivity. 
Assuming a decoupling of the rotational and translational motions they showed that in narrow tubes one can expect 
an increase of the selectivity in narrow tubes of a factor 100 at 20 K, when compared with models that approximate 
hydrogen as a sphere. 

These results have been confirmed by Trasca et al&, who calculated the D2 /H2 selectivity in the interstitial channels 
and groove sites of various carbon nanotubes bundles. 

Lu et al£ have also calculated, in the framework of the simple theory, the energy levels and selectivity of molecular 
hydrogen in carbon nanotubes, by numerical diagonalization of the Schrodinger equation of a confined rotor, finding 
a value of the total selectivity for T2/H2 mixtures at zero pressure and 20 K of the order of 100 in the (3,6) tube, 
a result far lower than the ones already published in the literature. This result was later attributed to the use of 
an unphysical potential for the hydrogen-carbon interaction. Calculations with more realistic potential s 10 ' 11 have 
confirmed the expectation of high selectivity. 

Two problems seem still without a definite solution: the first is the effectiveness of model potentials in describing 
the carbon- hydrogen interactions. Different authors have used different potentials, and their predicted values of the 
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selectivity scatter in an enormous range, from 10 2 to 10 7 . The second point, quite related to the first, is the actual 
effect of the rotational degrees of freedom on the selectivity. 

We have addressed both these issues in a recent paper—, where, in the framework of the simple theory, we have 
calculated the dependence of the zero pressure selectivity on the hydrogen-carbon interaction potential and we have 
developed an approximate method to evaluate the contribution of rotational-translational coupling to the selectivity 
itself. 

We have demonstrated that the zero pressure selectivity is very much dependent on even small changes in the 
interaction potential, especially in the case of narrow tubes. Moreover the presence of a steep confining potential results 
in a strong translational-rotational coupling, and the overall selectivity cannot be calculated neither by assuming 
independence between those two degrees of freedom nor by assuming a spherically symmetric model for the hydrogen 
molecule. These results have been validated by the exact diagonalization of the hydrogen-carbon Hamiltonian reported 
in a recent paper by Lu et al^ 

In this paper we address the issue of the contribution of rotational degrees of freedom in more detail, using the path 
integral Monte Carlo method to treat the quantized rotational degrees of freedom in an exact way. After discussing 
how to perform an efficient path integral simulation of rotors using the hybrid Monte Carlo (HMC) technique, we 
develop a method for the classical treatment of the rotational degrees of freedom only and show that it gives a much 
lower zero pressure selectivity than the full quantum model. We also extend the calculations at finite pressures and 
discuss how the selectivity changes in this regime. 

This paper is organized as follows. In Sec. |TT] we describe the hybrid Monte Carlo method that we have used in 
this work. We further address the issue on how the zero pressure selectivity can be calculated from path integral 
simulations, with particular emphasis to the influence of the quantized rotational degrees of freedom. We develop a 
formalism that enables one to simulate classical rotations together with quantized center of mass degrees of freedom 
using a path integral approach and use this method to investigate the effect of the quantized rotational degrees of 
freedom on the selectivity. 

In Sec. IIIII we present and discuss our results. Some technical derivations are presented in the appendices. 



II. THE HYBRID PATH INTEGRAL MONTE CARLO METHOD 
A. The path integral formulation of statistical mechanics 

The quantum mechanical expression for the partition function of a system of N rigid linear rotors is 

Q = Jd m X x d 2N n x (XiOxlexpI-^IXiOi) (1) 

where we denote with X a vector with the 3N center of mass coordinates of the rotors, and with O the set of the 2N 
angles describing their orientations. A subscript 1 has been introduced for later convenience. The Hamiltonian H is 
given by 

N h 2 N L 2 N 
H = X ~ 2M V * + X + X XjtJj) + Vext&Wi) (2) 

i—1 i—1 i<j i—1 

and is a function of the mass M, the angular momentum L and the inertia moment I of the molecules. We have 
introduced a fluid- fluid pair interaction potential v and a solid-fluid interaction potential w ox t- The quantities Xj 
and uji denote the operators for the center of mass position and the orientation of molecule i, respectively. The 
Hamiltonian is, as usual, the sum of the transnational (center of mass) kinetic energy T, the rotational kinetic energy 
K, the fluid-fluid potential energy V and the solid-fluid potential energy V eX t- 

The quantum partition function of Eq. ([T]) can be rewritten using repeatedly the Trotter identity 

exp[A + B}= lim (exp[A/P) exp[B/P]f (3) 

P— >oo 

valid for generally non-commuting operators A and B and further approximated by assuming a large but finite Trotter 
number P, obtaining 

r p 

Q~ / d m X 1 d 2N Q 1 (X 1 il 1 \l[exp(-pf/P)exp(-l3k/P)exp[-f3(V + V cx t)/P}\X 1 n 1 ) (4) 
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We can now introduce P — 1 completeness relations of the form 

1 = Jd 3N X l d 2N Q l \XitU)(XiSk\ (5) 
between the factors in Eq. (H|) and write the partition function as 

d 3N x 1 d 2N n 1 ...d m x P d 2N n P 



{XtSl^expi-pT/P) exp{-0K/P) exp[-0(V + V^ t ) / P]\X 2 tt 2 ) 
(X 2 n 2 1 exp(-0f/P) exp(-0K/P) exp[-0(V + V^/P] \X 3 n 3 ) 

(X P tlp\eM-0T/P) exp(-0K/P) exp[-0(V + V cxt ) / P^X^) 
d 3N X 1 d 2N tt 1 . . . d 3N x P d 2N n P 



p 

Y[(X t fl t | exp(-/3f /P) exp(-0K/P) exp[-0(V + V^t)/P] \X t+1 n t+1 ) (6) 
t=i 

where we have denoted Xp + \ = X\ and Qp+i = f2i. Each of the matrix elements appearing in the previous equation 
can be written as 

(XiSU\ exp(-/3f /P) eM~0K/P) exp[-/3(1/ + V ext )/P] \X i+1 Q i+1 ) = 
(X i \exp(-0t/P)\X i+1 ) exv(-0K/P)\il i+1 ) ex V [-0(V(X l+1 {l l+1 ) + V cxt (X l+1 n t+1 ))/P] (7) 

A straightforward calculation shows that the expectation value of the translational kinetic energy Boltzmann factor 
assumes the formic 

(Xi\exp(-0f /P)\Xi+i) = aexp(-/? K |X 4 - X l+1 \ 2 /2) (8) 
where the amplitude a and the "spring constant" k are given by 

(Mk B TP\ 3 / 2 (g) 

(!0) 



V 2Tth 2 

MP(k B T) 2 



h 

and the expectation value of the rotational kinetic energy Boltzmann factor becomes^ 

(SlileM-PK/PWi+i) = J2T,^^ P j( cos9 h+i) e M-W +VB/P] 

n=l J=0 
N 

= en) 

n=l 

where Pj(-) is a Legendre polynomial, is the angle between the directions Wj and LOi+i relative to molecule n 

and B = h 2 /(2I) is the rotational constant of the rotor. In the case of homonuclear molecules the indistinguishability 
of the nuclei imposes some restrictions on the sum in Eq. according to the spin states of the nuclei: for para-H 2 , 
ortho-D 2 and para-T 2 , as in this work, the summation on the angular momenta J in Eq. (|lip is limited to the even 
numbers onlyi^ and results in a positive definite density matrix, which can be directly used in the Monte Carlo 
simulations. 

In the other rotational states (i.e. ortho-H 2 and T2 and para-D 2 ) the sum is restricted to the odd angular momentum 
states, resulting in a density matrix which is not positive definite. As a consequence, more care has to be taken in 
performing a Monte Carlo simulation in this case^. 

The net effect of these algebraic manipulations is that we have been able to rewrite the original quantum partition 
function of an TV particle system, Eq. ([T]), as a classical partition function of a NP particle systems. The NP particle 
of the classical equivalent are naturally divided in P subsets (also known as time slices) of N particles each. Each 
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particle in a time slice interacts with all the other particles in the same time slice via the original intermolecular and 
intramolecular potential divided by a factor of P (see Eq. (J7J). Quantum mechanical effects taken into account by the 
interaction of each particle with the corresponding copy on the previous and following time slice: the center of mass 
coordinates are bound by the harmonic potential of Eq. © and the orientations give rise to the inter-slice rotational 
partition function of Eq. (jlip . The resulting system is then equivalent to a classical collection of N ring polymers, 
each having P beads. The i-th particle on a given polymer interacts only with the corresponding particle on the other 
polymers via the original intermolecular potential (rescaled by a factor of P) . The interaction between the beads of a 
given ring polymer are described by an harmonic interaction on the translational coordinates with the two adjacent 
beads (see Eq. (|8])) and an interaction between the orientational degrees of freedom of adjacent beads whose density 
matrix is given by Eq. 



B. Observables and estimators 



In the framework of the path integral approximation to the quantum partition function the estimators for the 
average values of interest can usually be obtained by their thermodynamic definitions. 
The estimator for the translational kinetic energy is theni^ 

T cst = \k B TP - f/ quant (12) 

where C/ qU ant = ( K /2)X)*Li l-^t — ^t+i| 2 is the value of the quantum spring potential energy, and the one for the 
rotational kinetic energy isi^ 



i p 

RC st = -J2K! st (13) 



i=l 



Kf st = ^E J(J+1)&J+1) Pj(coa9 iW )exp(-^BJ(J+l)j (14) 

We can derive an expression for the calculation of the particle density in a path integral simulation. Denoting x\ 
as the position operator of particle i in time slice t, one has 

N 1 

71— 1 

1 r N i 

= -J d ^X ld ^ (X 1 n 1 \Y / ^(r^xiP)exp[-pH]\X 1 n 1 ) 



(16) 



p 

(XrOxl f[ exp(-|f) exp(-|^) exp[-|(t> + V^X^) (17) 

introducing now P — 1 completeness relations, analogously to the passage leading to Eq. ([6]) and noting that the 
integrand can be rewritten by relabeling X\ <-> Xj, one obtains from Eq. p7p 



p / p N 



f[(X p n p \exp(-^f)eM-^K)eM-^(V + V cxt )]\X p+1 n p+1 ) (18) 
P =i 

so that the density at a point r is equivalent to the probability of finding the center of mass of a bead at the same 
point. 
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C. The hybrid Monte Carlo method 



Using the path integral formulation, a quantum partition function can be rewritten as a classical partition function 
of a system with a larger number of particles, so that classical Monte Carlo methods can then be used to calculate 
thermodynamic properties. Since we expect to work at conditions where quantum mechanics is not a small correction, 
it will be necessary to use large values of the Trotter number P. A simple Metropolis method for the sampling of the 
translational and rotational phase space, such as the one discussed in Ref. [14], will be affected by slow convergence. We 
have then decided to use the hybrid Monte Carlo method , 15 ' 16 that consists in choosing a new candidate configuration 
by performing a molecular dynamics (MD) move with a large timestep; the resulting configuration is then accepted 
or rejected using a standard Metropolis condition on the difference in the total energy (which is not conserved if a 
large enough timestep is taken). 

In order to perform an MD move one has to know the forces and the torques acting on each of the rotors. The 
potential energy between the molecules on the same slice is given by the rescaled original potential, and the quantum 
mechanical effects on the translational degrees of freedom are described by a simple harmonic potential between 
adjacent slices (see Eq. ([5])). The only unknown is the quantum torque between the molecules in adjacent slices: 
we have decided to calculate it numerically, starting from the expression of the density matrix in Eq. (jlip (see 
Appendix [B] for an analytic limit in the case of heteronuclear molecules) . Since Eq. (TlTj) represents a Boltzmann 
factor it can generally be written as B TOt — Cexp[— (3U ro t{9)] where C is an unknown constant and U ro t(9) is the 
quantum rotational potential energy between two adjacent rotors, whose orientations form an angle 9 with one another; 
we show in the appendix that for heteronuclear molecules in the high P limit U m t is given to a good approximation by 
the harmonic expression U ro t(9) = K9 2 /2 with K given in Eq. (|B3|) . In the general case the modulus of the quantum 
torque can be written as 

Equant [9) = — = ~^=- Sm6>- (19) 

d9 a a cos 9 

The direction of the torque is obviously orthogonal to the plane generated by the two orientations of the interacting 
molecules, and such that it tends to close the angle between the two molecules. The derivative in Eq. l[T9")) could be 
evaluated either numerically or using the identity 

dPijx) = IxPjjx) - IPj^jx) 
dx x 2 — 1 



D. The multiple time step method 

Inspection of Eqs. |0 and (jlOjl shows that the intermolecular forces scale like the inverse of the Trotter number P, 
whereas the quantum forces (and, possibly, the torques) are proportional to it. In order to efficiently sample phase 
space using an hybrid Monte Carlo method it is necessary that all the degrees of freedom contribute uniformly to the 
non-conservation of energy when a MD time step is performed. Since the intermolecular forces become weaker for 
large Trotter number while the quantum forces become stronger, we have decided to use a multiple time step method 
to perform the MD evolution^ We divide the forces into "long range" (the intermolecular forces, in our case) and 
"short range" (the quantum forces and torques) . Each of the long range dynamical steps At is divided into n smaller 
time steps St where only the short range forces are evaluated as the system is propagated. 

In our case we do not know the typical time scale of the quantum rotation, so we have decided to use three nested 
loops: we use a long time step At to propagate the system according to the intermolecular forces, an intermediate 
time step 5t to propagate the quantum spring forces on the translational degrees of freedom and a rotational time 
step St to propagate the quantum torques. 

We also need, for the hybrid method to work, a reversible algorithm to integrate all the degrees of freedom. It is 
well known that the velocity form of the Verlet algorithm possess such a feature and can be used in multiple time 
step methods. 17 Instead of using algorithms already developed to treat the general motion of rigid rotors in a multiple 
time step framework^ we have developed a velocity Verlet like integrator for the rotational motion of a rigid linear 
rotor that can be easily integrated in the velocity Verlet evolution of the center of mass coordinates. Details of its 
derivation are given in Appendix [Al 

Denoting by x and v the translational positions and velocities, and by e and w the direction of the molecular axes 
and the molecular angular velocities, the multiple time step method is then implemented as follows: 

Calculate Fi ong and Ni ong (intermolecular forces) 
v -»■ v + AtF long /(2M) 
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m^m + AtN long /{2I) 

Calculate F s f lort (quantum spring) 
v ^v + StF short /(2M) 

Calculate N s h or t (quantum torque) 
w — » w + 5TNi ong /(2I) 
e^e + Srmxe — (St) 2 vj 2 /2 
normalize e 

Calculate N s h ort (quantum torque) 
w -> w + 5TNi ong /(2I) 
x — » x + «(5r 

Calculate F s h ort 

v ^v + 5tF short /(2M) 
Calculate Fi ong and iVjong 
« -► « + AtF long /(2M) 
zu^zu + AtNi ong /(2I) 

Calculate the final translational and rotational kinetic energy 



In order to fix the integer values of the ratios 8t/8r and At/8t we have proceeded as follows. For the first ratio 
we have performed some ideal gas simulations, where the translational and rotational motions are decoupled and 
can be simulated by taking into account translations and rotations separately. We have adjusted the rotational 
and translational ideal gas time steps in order to have a 50% acceptance ratio and we have then fixed St /St — 

(<5t/5r) ideal gas — 3. 

The ratio At/8t is more difficult to set, but as a first guess one can set it equal to the inverse ratio of the Einstein 
frequencies corresponding to the intramolecular and quantum forces, that we calculate during the course of the 
simulations. For particles confined in narrow tubes we found that the optimal ratio At /St is of the order of 8-10. 
There are at least two possible strategies for the parallelization of the multiple time step algorithm for rigid rotors. 
In the first case, one can distribute the number of beads among the available nodes, and in the second case one 
can distribute the number of molecules. Using the first method requires very frequent and very short messages to 
be passed between the nodes, i.e. every short time step St. In the second case the quantum dynamics is performed 
locally, but the positions of all the molecules have to be redistributed among the nodes at every long time step in order 
to calculate the intermolecular forces: this results in less frequent communications (i.e. every long time step At), but 
the amount of data for a single communication is larger than in the first case. In the case of non-interacting molecules 
(i.e. the zero pressure limit) the last strategy results, of course, in a negligible amount of inter-node communication 
and is therefore optimal. 



E. Calculation of the selectivity 



The selectivity of two components, say T2 and H2 is defined as: 

XT 2 I XH 2 



S(T 2 /H 2 ) 



VT 2 /yH 2 ' 



(20) 



where x and y are the mole fractions in the adsorbed and bulk phases, respectively. 

In the limit of zero pressure, when one can neglect the adsorbate-adsorbate interaction, the selectivity depends 
only on the energy levels E of the adsorbed molecules, and can be written as^ 



S (T 2 /H 2 ) 



Qh 2 ° 



}%T Qh 2 



m H2 

VflTI 



d/2 ^jfree— rot 
Wt 2 



o tr< 



free — rot 



^exp(-^ 2 /fcT) 
_i 



(21) 



where <5 frec is the molecular partition function for the ideal gas, Q frec ^ rot is the free rotor molecular partition function 
and Q is the molecular partition function for the given specie. We have denoted by d the number of spatial dimensions 
in which confinement takes place. In the case of hydrogen molecules in carbon nanotubes, d = 2. In the zero 
pressure limit the selectivity is a function of the energy levels of the two species, which can be obtained by a direct 
diagonalization of the single-particle Hamiltonian. We shall refer to this procedure as "the simple model" . 
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We notice that, under the assumption that the rotational and translational dynamics are independent, one can 
rewrite Eq. (|2"Tj) as 

/ /ntras \ / /"}free— rot ^jrot \ 

C (V IVf \ ~ f m ^2 ^T 2 \ / ^T 2 VT 2 \ ctrasorot / ooA 

5o(T2/ff2) -l^Q^J^Q=QfJ- 5 ° *° (22) 

The calculation of the selectivity of a confined T2 /H 2 can in principle be performed using the definition given in 
Eq. (|2H)) . evaluating the mole fractions in the bulk and in the confined system using Grand Canonical Monte Carlo 
simulations. 

Since we expect quite high selectivities (possibly of the order of 10 5 or more), the evaluation of the mole fraction 
would require the simulation of very large systems, of the order of at least 10 5 molecules. In order to avoid the use 
of such demanding calculations, we have employed the method developed by Challa et al&^ to evaluate the zero 
pressure selectivity. With a straightforward extension to rigid rotors the selectivity in the limit of zero pressure can 
be written as 

So(T 2 /H 2 ) = C (exph/?A[/ H2 ^T 2 ])H 2 (23) 

where 

At/ H ^T 2 = r 2 dm (^) (24) 



IIIH'l 



V dm 



is the variation of the quantum potential energy when a H2 molecule is gradually transformed into a T 2 molecule by 
performing a number Nmc of Monte Carlo steps and the constant C is given by 



)frcc-rot( T2 ) 



where we have denoted by / the inertia moment of a given specie, and the average in Eq. ((23 j) is performed on a 
simulation of the lightest specie only. A number of the order of Nmc — 4000 points are necessary to reach convergence 
in the evaluation of the integral in Eq. (|24|) for the H 2 to T 2 transformation. 

In order to calculate the selectivity at finite pressures Challa et al~ have developed an efficient method to perform 
Grand Canonical simulations of mixtures. The usual insertion and deletion moves are performed for the heavier specie 
only (T 2 in our case), whereas the lighter specie is inserted or deleted by performing T 2 <-* H 2 transformations. A 
transformation move of a molecule of the specie 1 into a molecule of the specie 2 is accepted with the probability 



■Pi—,2 = mm 



1, ( y ) ex P^ 2 " ^ exp[-/3Af/ cxt ] 



(26) 



where N is the number of molecules already present in the system, A is the de Broglie wavelength, /j, the chemical 
potential (which is different for the two species to account for a given bulk molar composition) and AU ext is the 
difference of the sum of the fluid- fluid and solid- fluid potential energies between the configurations (JVi, A 2 ) and (Ni — 
1, N% + 1). Note that, in order to fulfill the detailed balance condition, the probability of a T 2 — ► H 2 transformation 
attempt must be equal to the probability of attempting the reverse move. 



F. Classical treatment of the rotational degrees of freedom 

In order to assess the importance of the quantized rotational degrees of freedom in quantum sieving we now develop 
a formalism to describe a system in which only the translational degrees of freedom are quantized and the rotations 
are described classicaly. 

In what follows we perform the derivation referring to a single rotor in an external potential, in order to avoid a 
cumbersome notation. The extension to interacting rotors is straightforward. 

Consider a system whose Hamiltonian is given by H = T(p) + K(L) + V(x, ft), where T(p) = p 2 /2m is the kinetic 
energy of translation, K(L) = BL 2 is the kinetic energy of rotation and V(x, O) is the potential energy with a 
dependence on the position operator x of the center of mass and the direction O. The quantum mechanical partition 
function is given by 



I. , m 



(27) 
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The classical treatment of some of the degrees of freedom correspond to the assumption that the operators of 
the corresponding generalized coordinates and momenta commute. Since we are interested in approximating the 
rotation as classical, we proceed as if the rotational kinetic energy and the potential, which depends on the molecular 
orientation, obey the commutation relation [V, K] = 0, which implies exp[— (3(T+V+K)] = exp[— /3(T+V)] exp[— f)K\. 
One can then perform the partial trace over the rotational degrees of freedom, obtaining 

q = Tr rot exp[-/3(T + K + V)] 

= ^2(l,m\exp[-(3(T + K + V)}\l,m) 

Im 

= ^2(l,m\exp[-0(T + V)]\l,m) e~ pm ^ 

Im 

= J dujid^Y^ (Jm|wi)<wi|e3cp[-)9(T + V0]|«2)<W2|J,m> e ~ f)Bl{l+1) 

Irn 

= J dwi<hj 2 '^2 ^m( w i)^m(^2)(wi|w 2 ) exp[-/?(T + V(x, uj 2 ) 



e -0Bl(l+l) 



-/*»(? 



I 

21 + 1 



4tt 



e -0Bi(i+i) exp [-j3(T + V (x,uj 1 )] 



f d 

= Qrot J -^exp[-p(T(p)+V(&,u))] (28) 

where Q ro t is the molecular partition function of the free rotor. In the last expression the direction to is the classical 
direction of the rotor and p and x are the momentum and position operators (still quantum mechanical). 

When one applies the Trotter formula to the reduced density matrix q each of the beads corresponding to a given 
molecule has the molecular axis pointing in the same direction as the others, since to is now considered as a classical 
variable. 

III. RESULTS AND DISCUSSION 
A. The potential model 

In order to assess the importance of quantized rotation on the selectivity, we need a potential model that explicitely 
treats the hydrogen molecule as a rigid rotor. To the best of our knowledge no such model has been extensively tested 
in the literature. We have recently evaluated the zero pressure selectivity of various potential models in the (3,6) 
nanotube in the framework of the "simple theory" , that is, using Eq. (|2"Tj) with the energy levels obtained by a direct 
diagonalization of the Hamiltonianj 10 ! 11 Many models predict very high selectivities (in the range 10 7 — 10 10 ), and we 
noticed a strong dependence of the selectivity on the potential parameters. 

We have chosen to use a reasonable potential model, well aware of the fact that it might not be the best one to 
describe the actual interaction of hydrogen molecules with themselves and with a carbon nanotube. We would like to 
point out that our main interest in this work is to evaluate the effect of quantized degrees of freedom on the selectivity, 
and not to develop an accurate model to describe the hydrogen-carbon or hydrogen-hydrogen interaction. 

We describe the hydrogen molecule as a rigid rotor of length / = 0.74 A with two Lennard- Jones sites on the position 
of the hydrogen atoms, having as parameters e = 8.4 K and a — 2.81 AiSi The Lennard-Jones radius of this model is 
less than the Lennard-Jones radius of the Buch potential (where CTBuch = 2.96 A), which has been shown to be a good 
spherical potential to describe the bulk properties of hydrogen. The smaller a takes into account the fact that the 
potential has two centers, separated by a distance corresponding to the actual gas phase H2 bond length. The value 
of e is almost one fourth of the corresponding value for the Buch potential (where eeuch = 34.2 K) and this takes into 
account the fact that in a site-site model we have actually four interactions between the two molecules. 

We have also described the carbon atoms in the nanotubes with a Lennard-Jones potential, using the Steele 
parameters ac = 3.4 A and ec = 28.0 Ki2i Solid-fluid interactions have been calculated using the Lorentz-Berthelot 
mixing rules. We have generated carbon nanotubes of various sizes and tabulated the solid- fluid potential by averaging, 
in cylindrical coordinates, over the length of a unit cell in the direction z of the tube axis and over the angle <f> in a 
plane orthogonal to the tube axis, thus obtaining the solid-fluid potential as a function of the distance of the molecule's 
site from the nanotube axis only. 
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12 3 4 

Radial distance [A] 

FIG. 1: The average potential energy curves for one of the hydrogen sites with the carbon nanotubes used in this study, as a 
function of the distance from the tube axis. 
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17 


-1313 


98 




(6,6) 


2.3 


-1016 


66 





(6,6) classical 


2.3 


-1018 


55.7 




(6,6) simple model 


2.12 


-1022 


41.9 




(10,10) 


5.6 


-524 


71 





(10,10) classical 


4.2 


-532 


66 




(10,10) simple model 


4.7 


-556 


58.1 





TABLE I: The values of the para-T2 / para-H2 selectivities calculated for different nanotubes at 20 K. The potential and kinetic 
energies reported are those of the H2 molecule. The simple model results have been obtained by assuming that a molecule is 
aligned with the nanotube axis and that rotations do not contribute to the selectivity. 

In this study we have focused our attention on the carbon nanotubes called (3,6), (2,8), (6,6) and (10,10) with the 
standard nomencalture. These tubes have geometrical radii of 3.1, 2.6, 4.1 and 5.1 A, respectively. We report the 
profile of the potential energy with one of the hydrogen sites if Fig. [1] 

B. The zero pressure selectivity 

Our results for the influence of the quantization of the rotational degrees of freedom on the selectivity are shown 
in Table HI where we report zero pressure selectivities and average energies for hydrogen confined in the (3,6), (2,8), 
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(6,6) and (10,10) carbon nanotube. 

The simulations have been performed using the method described in Eq. (|23[) . A system of at least 50 non 
interacting molecules was equilibrated inside a carbon nanotube for at least 20000 HMC steps, and the selectivity 
was then calculated by performing a production run of at least 20000 HMC moves. The selectivity was calculated 
from Eq. (|23|) , with configurations sampled every 50 HMC steps and using Nmc = 4000 points for the integration of 
Eq. (pi. 

It is apparent that the explicit inclusion of the rotational degrees of freedom has a dramatic effect on the selectivity, 
especially in the narrowest tube, where it jumps from 27000 up to 6.34 x 10 7 , a more than 2000 fold increase. In the 
other tubes, the effect of the rotational degrees of freedom is less dramatic, being of the order of 2.5 in the (2,8) tube 
and almost negligible in the (6,6). 

This enhancement is more than can be expected from a simple analysis that assumes independence of the trans- 
lational and rotational degrees of freedom^ 7 - which would give, in the case of the (3,6) tube with our potential, a 
selectivity 750 times higher than the one obtained assuming an isotropic potential. 1 - Rotational-translational cou- 
pling effects must be taken into account in narrow carbon nanotubes. We would also like to point out that the actual 
magnitude of the rotational-translational coupling on the selectivity calculated with the path integral simulations is 
still higher than the one obtained using the approximate method that we have developed in Ref. [lfj| . which would 
give, for the same system analyzed here, an enhancement of the selectivity by a factor of 1650. 

The physical origin of the higher selectivity due to the quantization of the rotational degrees of freedom can be 
seen by analizying the difference between the simulations in which the rotational are treated classically and the ones 
in which the rotations are treated quantized. 

In the presence of narrow confining potentials one expects to find the molecules aligned with the nanotube axis, 
when rotations are treated classically. In fact, we have calculated in Ref. [h] the selectivity in the framework of the 
simple model assuming a perfect alignment - so that the potential energy is a function of the center of mass position 
only - and we can see, from the results reported in Table HI that the results are in very good agreement with the path 
integral simulations where rotational degrees of freedom are assumed to behave classically. 

We have also calculated the average angle 9 between the molecular axis and the nanotube axis in the classical 
simulations as a function of the distance of the molecule from the nanotube axis. One can see from Fig. [5] that in the 
classical case one obtains O ~ 8 degrees, thus confirming an almost perfect alignement. 

The situation changes dramatically upon quantization of the rotational degrees of freedom, since the confining 
effect of the potential energy is now counterbalanced by the quantum delocalization. The average angle between 
the "orientation" of a given bead and the nanotube axis averages around 35 degrees at the center of the tube and has 
a slight decrease towards 25 degrees for the molecules which happen to be off center, where the potential energy is 
steeper. 

The effect of quantum fluctuations in the orientation can also be seen by the high value of the rotational kinetic 
energy reported in Table |U Hydrogen molecules confined in the narrowest tube have an average rotational kinetic 
energy up to 109 K, a very high value when compared to the free-rotor rotational energy at the same temperature, 
which is about 0.1 K due to the freezing of the rotational degrees of freedom on the I = spherically symmetric 
ground state. Due to the confining potential that tends to localize the molecular direction along the nanotube axis, 
the rotational state of the adsorbed rotor is a superposition of higher angular momentum states, resulting in a non-zero 
value of the average kinetic energy. 

If we now consider, in a semiclassical picture, a molecule at a given distance from the nanotube axis, we see that 
the quantization of the rotational degrees of freedom (and the consequent rotational delocalization) has the effect 
that the molecule samples region where the potential energy is higher with respect to an almost perfectly aligned 
(classical) configuration. One can then see that the average potential acting on the center of mass is steeper when 
rotations are quantized than in the case when rotations are treated classically, and that the steepness is higher for 
the lighter molecule than for the heavier one. As a consequence, the energy levels of the lighter rotors are more 
separated with respect to the energy levels of the heavier specie, not only because of a different mass, but also because 
the quantization of rotations has a different effect on the two kind of molecules. In the light of the simple theory 
formula, Eq. (|2ip , we can then see that quantized rotations enhance the spacing between the energy levels with a 
greater effect on the light isotope, thus having a big effect on the selectivity, as is indeed observed in the simulations. 
This phaenomenon has been termed "extreme two dimensional confinement" by Lu et al<^ By a careful analysis of 
numerically exact eigenstates of hydrogen isotopes confined in carbon nanotubes, they found that a large value of the 
selectivity is to be expected in geometries so narrow that the rotational ground state takes contribution from states 
with finite angular momentum. We have been able to show^ using an approximate model for the description of the 
coupled rotational and translational degrees of freedom, that under these circumstances a very large contribution to 
the selectivity does indeed come from the rotational degrees of freedom, as is apparent in the exact result that we 
show in Table U for the (3,6) tube. In the other tubes the average rotational kinetic energy has been found to be zero 
within the error bars, and the contribution of the rotational degrees of freedom to the selectivity is correspondingly 
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Distance from the tube axis [A] 

FIG. 2: The average angle between the nanotube and the molecular axes as a function of the distance from the (3,6) nanotubc 
center for H2. Solid line, full quantum simulations. Dashed line, classical treatment of the rotational degrees of freedom. 

smaller. 

One can also see the effect of the more confining effective potential by plotting the expected density of a molecule 
as a function of the distance from the tube center, when rotations are treated classically or quantized. We show the 
results for hydrogen in the (3,6) tube in Fig. [3j where it can be seen that the classical treatment of the rotational 
degrees of freedom results in a lower density along the tube axis, a signature of an effectively less confining external 
potential, as is also apparent from the average kinetic and potential energies reported in Table [TJ the quantized 
treatment of the rotational degrees of freedom results in higher average kinetic and potential energies, when compared 
to the case in which rotations are treated classically. 

C. Adsorption isotherms 

Using the path integral simulations we have calculated the adsorption isotherms for different hydrogen isotopes in 
the (3,6), (2,8) and (6,6) carbon nanotubes. The results are shown in Figs. [4] and [5] 

As a general trend, we notice that the heavier specie is the one most readily adsorbed, as is expected. Qualitatively, 
one can think that the larger thermal de Broglie wavelength of the lighter specie results in a larger effective Lennard- 
Jones radius, thus hindering the adsorption. 

Quantum effects result in a separation in pressure between the isotherms of different isotopes. This is a negligible 
effect in large tubes, such as the (6,6) where quantum mechanical effects do not influence the adsorption very much 
(as apparent from the low selectivity). In this case the isotherms of the different isotopes are very close, as can be 
seen in Fig. 0J 

The opposite is observed in the very confining (3,6) tube. It can be seen in Fig. [4] that T 2 can be adsorbed already 
at 10 -15 bar, whereas the adsorption of H 2 is not detectable below 10~ 5 bar, a difference of 10 orders of magnitude. 

In the (2,8) tube we observe a sizeable presence of quantum effects, though not as dramatic as in the (3,6), the 
isotherms of the lightest and heaviest species are separated by two order of magnitude in pressure. 

These isotherms allow us to clarify in detail how low the pressure must be in order to fulfill the zero pressure limit 
discussed in the previous section. For a given nanotube we can take as an estimate the pressure for which the isotherm 
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Distance from the tube axis [A] 

FIG. 3: Normalized densities for hydrogen confined in the (3,6) tube at 20 K from Eq. (|18jl . Solid line, quantum treatment of 
the rotational degrees of freedom. Dashed line, classical treatment of the rotational degrees of freedom 

of the heavier specie starts to be significantly different from zero. 

We notice that a high zero pressure selectivity is correlated with the separation in pressure of the isotherms 
corresponding to the two species. When the heavier molecules begins to be adsorbed, the lighter one cannot easily do 
so, hence the high values of the selectivity. 

D. Selectivity at finite pressures 

Using path integral simulations it is possible to investigate what happens to the selectivity at finite pressures. We 
have used the very efficient method of Eq. (f2"rl|) to calculate the selectivity using the definition in Eq. (|20[) . Differently 
to what happens in the zero pressure limit, the selectivity at finite pressure does depend on the assumed mole fraction 
in the bulk. We have assumed that at the low pressures of our study, the chemical potentials \i\ and fi2 of Eq. (|26|) 
can be evaluated using the ideal gas expressions. They are therefore a function of the total pressure and the assumed 
mole fraction in the bulk. 

Since we are calculating the selectivity using the observed mole fraction in the simulation box, it is apparent that 
we have to consider thermodynamic conditions where we expect finite amount of both isotopes to be present in the 
simulation box. For a given pressure we can then fix the bulk mole fraction so that the two species can be expected to 
be adsorbed more or less in the same amount. It is clear that for a given zero pressure selectivity Sq(A/B) bulk mole 
fractions such that ua/vb — l/S'o would result in an almost 1:1 ratio in the simulation cell. As a consequence less 
statistical errors can be expected when evaluating the ensemble average of Eq. (|2"0|) . We have decided to investigate 
the particular case of the (2,8) tube. Since observed a value of So(T2/rl2) = 52 (see Tabled, we fixed the bulk mole 
fraction of T 2 to the value yr 2 = 0.1. The results are shown in Fig. [6j 

We have also evaluated the selectivity in the framework of the ideal adsorbed solution theory (IAST), as outlined 
in Ref. 0. Assuming ideal solution behavior of the two species, it is possible to derive an expression of the selectivity 
starting from the isotherms reported in Fig. [5] We have fitted our isotherms with the same functional form used in 
Ref. H, and the quality of the fit can be seen on the same figure. The results of the IAST imply that the selectivity 
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Pressure [bar] 



FIG. 4: Pure fluid adsorption isotherms for hydrogen isotopes in the (3,6) (filled symbols) and (6,6) (open symbols) carbon 
nanotubes at 20 K. The circles and diamonds refer to para-Hb and para-T2 respectively, whereas the triangles refer to ortho-Y>2- 
The lines are drawn as a guide for the eye. 

raises monotonically with the pressure, and this prediction is related to the fact that upon saturation the density of 
the heaviest specie is slightly larger than the density of the heavier one. 

Also the path integral simulations show an increase of the selectivity at finite pressure, but this increase is not as 
steady as the IAST prediction. In fact the selectivity seems to tend to a constant value as the pressure is raised, as 
can be seen in Fig. [BJ We notice that the correct selectivity starts to deviate from the approximate IAST value at the 
pressure (10~ 21 bar in this case) where the adsorption isotherms begin to show saturation. 

The reason of this behaviour can be traced back to the structure of the adsorbed phase. Each of the two pure 
species occupy the same nanotube with an almost equal linear density around 0.27 molecules/ A. Therefore one might 
expect an "average distance" L ~ 3.7 A between the molecules in the case of a saturated mixture. 

In going from zero to finite pressures one can then expect that the motion along the z coordinate of a given molecule 
is progressively hindered by the ever close presence of the other ones, until a saturation condition is reached, and 
further compression of the system becomes more difficult. 

This picture can be validated by the results reported in Fig. [7J where we plot the average center of mass kinetic 
energies for H2 and T2 as well as the average fluid-fluid energy for the mixture adsorbed in the (2,8) tube. 

We notice that the average potential energy per molecule tends to a constant value of V ~ —16.5 K, indicating 
that the molecular configuration does not change very much after the onset of saturation. 

The effect of the localization along the z direction is apparent from the behaviour of the kinetic energies as a 
function of the loading. The values for T2 and H2, reported in Fig. do indeed show an increase at finite loading 
with respect to the zero pressure value. This indicates a progressive hindrance of the motion along the nanotube axis, 
until the saturation point is reached and the kinetic energy tends to a constant value. 

It is possibilc to make an estimate of these effects, by assuming that - at saturation - all the molecules in the tube 
are separated from their neighbors by the same distance, L m i n . In order to calculate L m i n one can assume that the 
molecules interact with an effective potential obtained by averaging the original potential on the rotational degrees of 
freedom. This takes into account the fact that the (2,8) tube is wide enough so that the molecules are freely rotating. 
It is then possible to calculate the average potential energy per particle, Vff{x), as a function of the distance x 
between the molecules, assuming that every molecule has two nearest neighbors located at ±x. 
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FIG. 6: Finite pressure T2/H2 selectivity in the (2,8) tube at 20 K, assuming a bulk mole fraction of T2 equal to yr 2 = 0.1. 
Solid line, prediction from IAST theory. Circles, path integral simulations. 
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FIG. 7: Average translational kinetic energy (circles, H2; triangles, T2) and average solid fluid potential energy (diamonds) for 
the adsorption of a T2/H2 mixture (j/t 2 = 0.1) in the (2,8) carbon nanotube at 20 K. 



The distance L min is then calculated as the point of mechanical equilibrium, 8Vff{x)/ dx = 0, obtaining the value 
Lmin — 3.43 A, not too distant from the experimental average distance given above. In this configuration the potential 
energy per particle is given by VFF(L m in) = —51.6 K. 

We further assume that any molecule performs harmonic motions around the potential energy minimum. In the 
actual system, of course, the dynamics is determined by anharmonic effects, so that the following estimate has only a 
heuristic value. The "spring constant" for these oscillation is evaluated to be kno = d 2 VFF/dx 2 (L min ) = 316.8 K/A 2 . 
Hydrogen isotopes oscillating in this harmonic potential would have a zero point energy of Eq 2 = tiy/kno /^Ti/^ = 
25.2 K and E^ 2 = 43.7 K for T 2 and H 2 respectively, so that we might expect V > V F F(L min ) + E„ 2 = -26.4 K, 
close to the actual value observed. 

One can further estimate the asymptotic value of the selectivity at finite pressure by assuming the that its value 
at saturation is given by the product of the zero pressure value Sq and the selectivity Sho due to the effective ID 
harmonic oscillator discussed above. Using and independent particle approach^ the value of Sho corresponding to 
T2 and H2 in a harmonic oscillator of spring constant kno turns out to be Sho = 1-6, in reasonable agreement with 
the observed behavior of the selectivity which tends, at saturation, to twice the zero pressure value. 



IV. CONCLUSIONS 



In this paper we have presented an algorithm for the hybrid path integral Monte Carlo simulation of rigid diatomic 
molecules. We show how to calculate torques from the expression of the rotational density matrix and we moreover 
show an approximate expression for the rotational density matrix of heteronuclear molecules that, in the large P 
limit, is completely analogous to the expression for the translational degrees of freedom. We have developed a velocity 
Verlet like integrator for the rotational degrees of freedom suitable for a multiple time step approach. 

We have then applied this method to the calculation of the para-T2 / para-Yi.2 selectivity in various carbon nanotubes 
at 20 K. We have discussed, in particular, the effect of the quantized rotational degrees of freedom on the selectivity by 
developing a simulation method that allows a classical treatement of the rotations while keeping a quantum tratment 
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of the translational degrees of freedom. 

We show that the explicit inclusion of quantized rotations enhances the zero pressure selectivity by a factor of more 
than 2000 in the narrowest tube we have investigated (the (3,6)), but does not have such a dramatic effect in larger 
tubes. The rotational degrees of freedom contribute by less than a factor of two in the (2,8) and larger tubes. 

We have been able to investigate the effect of finite pressures on the adsorption and the selectivity. We have 
calculated the adsorption isotherms of various hydrogen isotopes in different tubes at 20 K and shown that quantum 
effects hinder the adsorption of the lighter specie, whose isoterm can be separated by many order of magnitude in 
pressure from the one of the heavier specie. 

We have been able to calculate pressure dependence on the selectivity on the pressure and found that, for a 0.1 
molar mixture of T 2 in the (2,8) tube, the selectivity tends to twice its zero pressure value when the pressure is raised. 
We correlate this behaviour to the hindrance of the molecular motion along the nanotube axis. 
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APPENDIX A: VELOCITY- VERLET INTEGRATOR FOR RIGID LINEAR MOLECULES 



The dynamics of a rigid rotor is described by the equations 

de 



dt 
~dt 



ru x e 



= N 



(Al) 
(A2) 



where e is a unit vector in the direction of the rotor and, w is the angular velocity, I is the inertia moment and N is 
the torque applied to the system. 

These dynamical variable are redundant, since it can be seen that the norm of the vector e is a constant of the 
motion described by the previous equations. Since the torque N is, by construction, always orthogonal to the axis 
vector, the component of the angular velocity along the unit vector is also a constant of motion, and is usually set as 
zero. 

The previous equations cannot be then put in the form of a Hamiltonian system. In order to develop a time- 
reversible integrator that can be used in the hybrid Monte Carlo method, we demonstrate that it is indeed possibile 
to integrate the equations using a velocity Verlet like, adapted to take into account the abovementioned constraint. 
The multi-step algorithm can them be developed by analogy to the velocity Verlet case. 

Using the Taylor expansion we can write 



e(t + St) 



r de 1 , p n9 d 2 e 
^e(t) + St- + -(St) 2 - 

= e + St(zu x e) + ^(^ t ) 2 ^( w x e ) 



8t(w x e(i)) 



i(5t) 2 (~T' x e + nj x (w x e) 



e(i) + St 



w(t) 



St 



(sty 



■m 2 e(t) 



(A3) 
(A4) 
(A5) 

(A6) 



where we have restored the explicit time dependence on time in the last passage. 

The last term in Eq. (|A6|) assures that the length of the unit vector describing the direction of the rotor remains 
fixed. In our code we re-normalize the unit vector after each time step. 

The equation for the angular velocity becomes 



w(t + St) 



, . . dm 
w(t) + St— 



(St) 2 d 2 



dt 2 



St . St.. r .. 
= w + ~~^ w + ~~2 \ m ' Stvj) 



w(t) 



St N(t) St N(t + St) 



2 12 I 

So that one can construct a velocity Verlet like algorithm for the rotational degrees of freedom as 

1. Calculate the angular velocity at half time step w(t + St/2) — w(t) + ^-^7-^ 

2. Advance the orientation e at full time step, using Eq. (|A6|) 

3. Calculate the torques at the time t + St 

4. Advance the angular velocity at full time step 



(A7) 
(A8) 
(A9) 



APPENDIX B: APPROXIMATE TREATMENT OF THE ROTATIONAL DENSITY MATRIX 

It is known to be impossible to rewrite the rotational partition function in Eq. pip as the partition function of an 
effective harmonic potential, as happens for the translational degrees of freedom. 

We have noticed, though, that if P is large enough the rotational partition function of a hcteronuclear molecule 
can be approximated by an expression analogous to Eq. ([8]) for the translational degrees of freedom. 
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FIG. 8: The difference between the approximated and the numerical rotational density matrix for HD (solid line) and DT 
(dashed line) at 20 K with P = 100. The numerical density matrices for the same Trotter numbers are shown in the inset. 
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where 



A - -sr <B2) 

A* - ™^ (B3) 

Eq. (|B2p can be demonstrated by setting #y = and approximating the sum over J in Eq. (fTTj) with an integral. 
We have not been able to demonstrate Eq. (|B1[) analytically, but we have verified that the first 20 coefficients of the 
expansion in 9ij of both sides of Eq. (|Bip are the same in the P — > oc limit and that the exact and approximate 
density matrices are very similar in the same limit. We would like to point out the formal similarity between Eqs (|B2|) 
and (|B3[) and the corresponding ones for the translational degrees of freedom, Eqs. © and (fTT))) : in the large P limit 
the quantum mechanical effect corresponds to the action of an harmonic torque between the rotors in adjacent time 
slices. 

We show in Fig. [5] the difference between the numerical and the approximated rotational density matrix of Eq. (|B1[) . 
for HD and DT at a temperature T = 20 K, using a Trotter number P = 100. We can see that the approximate 
expression is accurate within 1.5% in the low region for the ligher specie. At a constant Trotter number, the 
approximation is of course better for a heavier molecule such as DT. 

One can see from Eq. (|Bip that in the large P limit the quantum rotational density matrix is analogous to the 
translational case (see Eq. (|Sjl), i.e. represents an effective harmonic potential energy between adjacent beads in the 
classical mapping. 

This is also apparent in the calculation of the torques, that we use in the hybrid Monte Carlo calculation. We 
report in Fig. [5] the results obtained using the numerical method derived in Eq. (JTHJ) and the analytical calculation, 
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FIG. 9: The quantum torque calculated numerically using Eq. (|19[) (symbols) and analytically using the "spring constant" 
given in Eq. (|B3|I (lines). Solid line and circles refer to HD at 20K, dashed line and triangles to DT at 20K. In both cases we 
have set P = 100. 



using the "spring constant" of Eq. (|B3[) . The two methods give virtually the same values, even though the numerical 
calculation shows numerical instabilities for angles just a little higher than the one shown in Fig. [5] 

Unfortunately the approximation of Eq. (|B 1 1) is not applicable for ortho- and para- species, because Sij(9) does not 
fall monotonously to zero for large 9, as implied by the right hand side of Eq. (|B1|) . 

If the summation in Eq. (|B1[) is restricted to the even angular momenta only, one obtains a density matrix with 
the propriety Sy-(7r/2 + 9) = By(7r/2 — 9), such as the one shown in Fig. [TUJ (Sy-(#) would be odd with respect to 
the point 9 = tt/2 if the odd angular momenta only were used in Eq. (|B1|) .) 

The approximation of Eq. (|B1[) is still valid for 9 < tt/2, with the modification A — > A/2. Since the density matrix 
is even with respect to 9 = tt/2, one could try to approximate it as 



9 a) - ^exp 



K9l 
2k B T 



@{tt/2 



A 

2 6XP 



K(tt - 9 K 



2k B T 



®(x-tt/2) 



where we have denoted by Q(x) the step function, which is for x < and 1 for x > 0. Unfortunately this will result in 
a discontinuous quantum torque at 9 = tt/2, and we have not pursued the idea further. Nonetheless the approximate 
expression of Eqs. (|B1[) - (|B3[) may prove to be useful for an efficient quantum simulation of heteronuclear molecules, 
where no such discontinuities appear. 

It is interesting to notice that in the classical limit Ti — > the torque constant of Eq. (|B3[) goes to infinity. One 
could then conclude, albeit heuristically, that a classical treatment of the rotational degrees of freedom only can be 
done by assuming that the orientations of all the beads corresponding to one molecule all point in the same direction, 
in agreement with the derivation presented in Sec. Ill Ft 

In order to test the validity of the approximation in Eq. (|B1[) . we have evaluated the zero pressure selectivity of a 
DT/HT mixture in the (3,6) tube using the approximate and exact method to evaluate the rotational density matrix. 
The results, shown in Table [IT] show that the approximate method compares very well with the exact one, even if 
the rotational kinetic energies are sensibly different (but still within the quite large error bars) and the zero pressure 
selectivity of a DT/HD mixture is overestimated by 50%. In this type of calculation, which considered N = 64 
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FIG. 10: Density matrix (solid line) and corresponding torque (dashed line) for para-H2 at 20K, using P = 100. 



Property 


Exact result 


Approximate result 


Solid-fluid potential energy (K) 


-745 ± 5 


-757 ± 5 


Translational kinetic energy (K) 


252.8 ±0.5 


253.2 ± 0.5 


Rotational kinetic energy (K) 


172 ± 43 


192 ± 44 


(Ad) (degrees) 


17.5 ±0.2 


17.4 ±0.2 


Rgyr (A) 


0.314 ±0.008 


0.316 ±0.008 


5o(DT/HD) 


5450 


8470 



TABLE II: Comparison of the approximate and exact method for the calculation of the density matrix. Properties of HD in 
a (3,6) carbon nanotube at 20 K with P = 100. (A9) is the average angle between the orientation of two adjacent beads and 
Rgyr is the gyration radius of the classical ring polymer. 



particles with P = 100 and parallclization on 8 nodes, we have observed a 35% speedup using the approximate 
method. 
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